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Abstract 

Segmentation plays an important role in many preprocessing stages in image processing. 
Recently, convex relaxation methods for image multi-labeling were proposed in the litera- 
ture. Often these models involve the total variation (TV) semi-norm as regularizing term. 
However, it is well-known that the TV ftmctional is not optimal for the segmentation 
of textured regions. In recent years directional representation systems were proposed to 
cope with curved singularities in images. In particular, curvelets and shearlets provide an 
optimally sparse approximation in the class of piecewise smooth functions with sin- 
gularity boundaries. In this paper, we demonstrate that the discrete shearlet transform 
is suited as regularizer for the segmentation of curved structures. Neither the shearlet 
nor the curvelet transform where used as regularizer in a segmentation model so far. To 
this end, we have implemented a translation invariant finite discrete shearlet transform 
based on the FFT. We describe how the shearlet transform can be incorporated within the 
multi-label segmentation model and show how to find a minimizer of the corresponding 
functional by applying an alternating direction method of multipliers. Here the Parseval 
frame property of our shearlets comes into play. We demonstrate by numerical examples 
that the shearlet regularized model can better segment curved textures than the TV reg- 
ularized one and that the method can also cope with regularizers obtained from non-local 
means. 

1 Introduction 

In recent years, much effort has been spent to design directional representation systems for 
images such as curvelets [8], ridgelets [9] and shearlets [28] and corresponding transforms 
(this list is not complete). Among these transforms, the shearlet transform stands out since 
it stems from a square-integrable group representation [13] and has the corresponding useful 
mathematical properties. Moreover, similarly as wavelets are related to Besov spaces via 
atomic decompositions, shearlets correspond to certain function spaces, the so-called shearlet 
coorbit spaces [l4j. In addition shearlets provide an optimally sparse approximation in the 
class of piecewise smooth functions with singularity curves, i.e., 

\\f-fN\\l,<CN-\logNf asiV^oo, 

where /at is the nonlinear shearlet approximation of a function / from this class obtained by 
taking the N largest shearlet coefficients in absolute value. 

In the following, we want to show how the directional information encoded by the shearlet 
transform can be used in image segmentation. Although shearlets and curvelets have been 
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applied to a wide field of image processing tasks, see, e.g., [22l [291 EQl EH [361 US], to the best 
of our knowledge they were not used in connection with segmentation so far. Segmentation 
is a fundamental task in image processing which plays a role in many preprocessing steps. 
Recently, convex relaxation methods for image multi-labeling were addressed by several au- 
thors [2l[l6l[26l[271[32l[33l[3ll[39l[ill[17]. In these methods a convex functional consisting of 
a data term and a regularization term is minimized. Due to its edge-preserving properties a 
frequent choice for the regularization term is the total variation (TV) semi-norm introduced 
by Rudin, Osher and Fatemi [ID]. The data term is either related to a predefined codebook (of 
gray/RGB-values related to the labels) or the codebook is updated during the minimization 
process itself. For the later approach see, e.g., [3l[TTl|25]. In this paper, we focus on the first 
approach, but replace the TV regularizer by the shearlet transform. To this end, we introduce 
a simple discrete shearlet transform which translates the shearlets over the full grid at each 
scale and for each direction. Using the FFT this transform can be still realized in a fast way. 
We show by numerical examples that the shearlet regularized model is well-suited to segment 
curved textures and can even cope with models incorporating nonlocal means. 



•1^ 



(a) Forms with different edge ori 
entations 




Figure 1: Shearlet coefficients can detect edges with different orientations. 

This paper is organized as follows: In Section [2] we introduce the translation invariant, finite, 
discrete shearlet transform whose full description cannot be found in the literature so far. 
Here we follow the path via the continuous shearlet transform, its counterpart on cones and 
finally its discretization on the full grid. This is different to other implementations as, e.g., 
in S'/ieariaiQ see [l3]. Our discrete shearlet transform can be efficiently computed by the 
fast Fourier transform (FFT). The discrete shearlets constitute a Parseval frame of the finite 
Euclidean space such that the inversion of the shearlet transform can be simply done by 



applying the adjoint transform. This is proven in Theorem 2.2 and cannot be found in the 
literature so far. In Section [3] we recall a convex multi-label segmentation model and show 
how the shearlet transform can be incorporated into this approach. We propose the method of 
alternating direction of multipliers (ADMM) Plll5 l [T7 1 ll8j to find aminimizer of the functional 
where the Parseval frame property of our shearlet frame can be exploited. Numerical results 
are presented in Section [4j in particular in comparison with models involving TV and nonlocal 
means (NL) regularizers [6l [19 ] [20 l lii] . 
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2 Finite Discrete Shearlet Transform 



We start with a brief introduction of the finite discrete shearlet transform to make the paper 
self-contained. First we describe the continuous shearlet transform in and on cones. Then 
we provide a discretization of the shearlet transform on a finite domain which is translation 
invariant and cannot be found in this way in the literature. We prove that the discrete 
shearlets constitute a Parseval frame of the finite Euclidean space so that the inversion of the 
shearlet transform can be simply done by applying its adjoint transform. For a more detailed 
description we refer to [23]. 



Continuous Shearlets. Let the parabolic scaling matrix Aa and the shear matrix Sg be 
defined by 

a \ ^1 „ /is 



^/a 



a G M+, Ss 



1 



s £ 



Then the continuous shearlets ^a,s,t emerge by dilation, shearing and translation of a function 
■0 e L2(M) as 

i;a,sA^) ■■= a-hiA-^S-\x - t)), (1) 
see |28l I31j . Using the Fourier transform T : L2(M^) — t- L2(M^) given by 

Tfiuj) = /(a;) := / /(t)e-2-<-'*)dt, 



we obtain that 

^a,s,t(w) = a3e-2'^*<^'*>V' (awi, ^/S(sa;i + ^2)) . 
The continuous shearlet transform SH^^f) of a function / € L2(M) is defined by 

S'Ui,{f){a,s,t) := {f,^a,s,t) = {f,1pa,s,t)- 
The shearlet transform is invertible if the function ip fulfills the admissibility property 

|V'(W1,W2)P 



dujiduj2 < 00. 



This is in particular the case if 



CO2 



where ^1 : R - 
functions are 



is a wavelet and ip2 '■ 



a bump function. Typical choices for these 



Tpiiioi) := y^b'^{2iOi) + b'^{uji) and -02(1^2) := 



a/u(T+7j^ for UJ2 < 0, 

y/v{l — UJ2) for UJ2 > 0, 



(2) 



where 



v{x) := < 



for X < 0, 
35x^ - 84x5 + 70x^ - 20x^ for < x < 1, 

1 for X > 1 
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and 



'sin(fv(|x| - 1)) for 1 < |x| < 2, 



b{x) := < 



^2 

cos(f7;(i|x| - 1)) for2<|x|<4, 







otherwise 



were introduced by Y. Meyer [381 [37] . The functions i/ji and ip2 are shown in Fig. [2j It is 
well-known that they have the following useful properties. 

Lemma 2.1. The functions tpi and Tp2 defined by ^ fulfill 

^|V;i(2-2^'u;i)|2 = l /or |u;i|>l, 



2^ 



\Mk + 2^W2)|' = 1 for \UJ2\ < 1, i > 0. 



k=-2i 




\ I 




Figure 2: The functions t/^i with supp^i = [—4, — |] U [^,4] (left) and 
ip2 with suppV'2 = [—1)1] (right). 



Shearlets on Cones. To get a good directional selectivity both from the horizontal and 
vertical point of view shearlets on the cone were introduced in ^28J. We define the horizontal 
and the vertical cone by 

■.= {{coi,co2) eR^ : Iwil > ^,\uj2\ < \u;i\}, 

C ■.= {{C01,C02) eR'^ : \UJ2\ > ^,\UJ2\ > \u}i\}, 

resp., and the "intersection" (seam lines) of the two cones and the "low frequency" set by 



1 



1 



:= {(wi, W2) G R^ ■■ |wi| > ^, |^^2| > ^, |wi| = |a;2|}, 
C° := {(a;i,W2) G R^ ■ |wi| < 1, |^^2| < 1}, 

resp., see Fig. [sj Altogether = U U U with an overlapping domain := 
(—1,1)^ \ {—^, 2)^- For each set C, k G {h,v, x} we define a characteristic function xc 
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Figure 3: The sets C^,r,C'' and C°. 



which is equal to 1 for G C and otherwise. The shearlets on the cone are defined by 



il^^{uji,uj2) ■■= V'i(wi)V'2 ( ^ j Xc^ and 'ip'" {ui,uj2) ■= '4^i{i^2)'>p2 



^^2 / 



Fig. [4] shows shearlets both in the time and frequency domain. Note that for |s| < 1 — -y/a we 
have supp '(/'a,s,t ^ and for 1 — ^/a < |s| < 1 + ^/a parts of supp-i/^^.s,* are also in C" which 
are cut off. For |s| > 1 + ^/a the whole shearlets are set to zero by the characteristic function. 
For {ijJi,0J2) G C^, i.e. \uji\ = \ijJ2\1 we define ■i/'^ (<^i, ^2) := V'('^i) '^2)xcx • Finally, we use the 
scaling function (j) defined by 



{^1,^2. 



{^p{oJi) for \uji\ < l,\u}2\ < |wi|, 
ip{u}2) for \UJ2\ < < \U}2\ 



with 



for < ^, 



cos(|u(2|(u| - 1)) for I <\u)\ < 1, 
otherwise 



and its translates (/>t(x) = 4'{x — t), i.e., 4>ti^) = e on Note that this scaling 

function (p (respectively if) matches perfectly with ipi, more precisely, since 



j>0 



sin^ (^v 



for < 5, 
2 .-(20; - 1)) for ^ < |cj| < 1, 
for \u}\ > 1 



we obtain 



\M^)\^ + \^H\^ = '^ for G 



(3) 
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Finite Discrete Shearlets. In the following, we consider digital images as functions sam- 
pled on the grid -^X := -^{(mi, 771,2) : = 0, . . . , — 1, z = 1,2}, where we assume 
quadratic images for simplicity. Further, let jo := [2log2A^J be the number of considered 
scales. To obtain a discrete shearlet transform, we discretize the scaling, shear and translation 
parameters as 

aj:=2'^^ = ^, j = 0,...,io-l, 



m 



-2^ <k<2\ 
m G I. 



With these notations our shearlets becomes tpj^k,mix) '■= '4'aj,Sj f,,tm,{x) = ip{A^^Sg_,^^{x — tm))- 



Observe that compared to the continuous shearlets defined in (|1|) we omit the factor a 4 . In 
the Fourier domain we obtain 



where O := {(a;i,(x>2) : = — , . . . , [y] — 1, i = 1,2} . By definition we have a < 1 and 



r 



J 



(a) Shearlet in Fourier domain 
1 

" 2 




for a, — ^ and s = 



(b) Same shearlet in time domain (zoomed) 



Figure 4: Shearlet in Fourier and time domain. 



|s| < 1. Therefore we see that we have a cut off due to the cone boundaries only for \k\ = 2^ 
where |s| = 1. For both cones we have for |s| = 1 two "half" shearlets with a gap at the seam 
line. None of the shearlets are defined on the seam line . To obtain "full" shearlets at the 
seam lines we "glue" the three parts together, thus, we define for |A:| =2^ a sum of shearlets 



i.hxv 



We define the discrete shearlet transform as 



Sn{f){K,j,k,m) 



if, (km) 



for K = 0, 
for K G {h, v}, 
for K = x,\k\ 



23 . 
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where j = 0, . . . , jo — 1, —2^ + l<k<2^ — 1 and m G X if not stated in another way. The 
shear let transform can be efficiently realized by applying the f f t2 and its inverse if f t2 which 
compute the following discrete Fourier transforms with 0{N^ log N) arithmetic operations: 

We have the Plancherel formula ^ 

{f,9) = j^if^a)- 

Thus, the discrete shearlet transform can be computed for k = h as follows: 

Sn{f){h,j,k,m) = {f,ij'i,^j = -^{f,i''lk,J (4) 

= ifft2(5j-fc) 

and with the corresponding modifications for k E {v, x,0}. 

In view of the inverse shearlet transform we prove that our discrete shearlets constitute a 
Parseval frame of the finite Euclidean space L2(I). Recall that for a Hilbert space Ti a 
sequence {uj : j G J} is a frame if and only if there exist constants < A < B < oo such 
that 

A\\ff<Y,\{f^^3)\'<B\\ff for all /e^. 

The frame is called tight if A = B and a Parseval frame if A = B = 1. For Parseval frames 
the reconstruction formula 

/ = E^/'^i)^i for ah /e 7^ 

holds true, see [12^ I37j . In the n-dimensional Euclidean space we can arrange the frame 
elements Uj,j = l,...,h>nas rows of a matrix U. Then we have indeed a frame if U has 
full rank and a Parseval frame if and only if V^U = In- Note that ULf^ = la is only true if 
the frame is an orthonormal basis. The Parseval frame transform and its inverse read 

{{f,u,))%, = Uf and f = U-'{{f,nj))%,. (5) 
By the following theorem our shearlets provide such a convenient system. 
Theorem 2.2. The discrete shearlet system 

{^lk,m, i^lk,m, ^f±2^,m : J = 0, . . . , Jo " 1 , -2-'' + 1 < A; < 2^' - 1 , m G X} U {0„ : m G X] 
provides a Parseval frame for L2(X). 
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Proof. We have to show that 

io-l 2J-1 io-1 

ii/ii'= E E E EK/'^i,^,Ji'+E E EK/'^irJi'+EK/'^'")!'-^- 

By Q we know that 

so that by Parseval's formula 

E \if^^km)\" = E = = . 

•max mdX 

Analogous computations can be done for the vertical cone, the seam-line part and the low-pass 
part, see [Ml- Summing up the different pieces we obtain 

, /io-l 2^-1 

^ = ]^ E E Ei^(4"''^i'4~'^^i+2-^'^2)iv>i,^2)p 

y j=0 fc=-23+lwe!^ 
jO-1 2J-1 

+ E E El^(4"''^2,4-^fcu;2 + 2-^'u;i)|2|/('^i,W2)|' 
i=o fc=-2^+i'<^eO 
io-i / 

+ E E E m-'^iA-'kuji+2-^uj2)\''\f{io,,u2)Wc^ 

+ ^ \i,{A-^U2,4-^kuJ2 + 2-^UJi)\^\f{iOi,C02)\\c^ 

+ J^|Vi(4-V,4-^fc^i + 2-^'^2)lV>i,^2)Pxcx ) 
+ J]|,^K,^2)|'|/(a;i,^2)n 
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which can be further simphfied by Lemma 2.1 



as 



weC^ i=0 A;=-2J ""^ 



EElfor|a;i|>l =1 

+ 1/(^1, cus)^^ 1^1 (4-^'a;2)|' ^ 1^^2(2^- + 
ijeC" i=o fc=-2J 



= 1 for |aJ2|>l =1 
+ E E l/('^l'^2)|'| ^{i-'0Jl,0) I^Xcx 

Finally we obtain by ([3| that 

^ = ^{ E l/(^i>'^2)|' + |/(wi,a;2)|'(|^i(u;i)|' + |^iMl' + |<A(wi,o;2) 

= ]^E 1/(^1' ^^2)1' = ii/f. 

and we are done. □ 

Having a Parseval frame the inverse shearlet transform reads 

io-i 2^-1 jo-l 

/= E E E E(/'^.-^,"^)^i,^,r.+E E E(/'V''rJ^?r..+E(/''^-)'^-- 

Ke{/i,i;} i=0 fc=-2J+l meX j=0 fc=±2J meX meX 

The actual computation of / from given coefficients c{hi, j, k,m) := {f-,4''jkrn) done in the 
Fourier domain. Due to the linearity of the Fourier transform we get, e.g., for the horizontal 
cone 

fHXc>^=Y. E E^^'^Jfc,™)^iAmM 
j=0 fe=_2J+l meX 

= E E E^(^'-^'''^'"')^~^"^^'^'"'^^(^~'''^l'^^'^'^l+2~^'^2). 
j=0 fc=-2J+l meX 

V ' 

fft2{c(/i,j,fc,-))(aJi,a;2) 

The inner sum can be interpreted as a two-dimensional discrete Fourier transform and can be 
computed via a FFT. Hence, / can be computed by simple multiplications of these Fourier- 
transformed shearlet coefficients with the dilated and sheared spectra of ip and afterwards 
summing over all scales j and all shears k. The spectra are the same as for the transform 
itself and can be reused. Finally we get / itself by an iFFT of /. For implementation details 
we refer to Ei 
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Remark 2.3. 1. We are aware of our larger oversampling factor in comparison with, e.g., 
ShearLab. Having 4 scales we obtain 61 images of the same size as the original image. But 
since shearlets are designed to detect edges in images we like to avoid any down-sampling 
and keep translation invariance. A possibility to reduce the memory usage would be to use 
the compact support of the shearlets in the frequency domain and only compute them on 
a "relevant" region. But we then would also have to store the position and size of each 
region what decreases the memory savings and would make the implementation a lot more 
complicated. 

2. In and the respective implementation ShearLab a pseudo-polar Fourier-transform is 
used to implement a discrete (or digital) shearlet transform. For the scale a and the shear s 
the same discretization is used. But for the translation t the authors set tj ^^^ ■= AajSg i^m. 
Thus, their discrete shearlet becomes 

Since the operation SJ. ^AajUJ would destroy the pseudo-polar grid a "slight" adjustment is 
made and the exponential term is replaced by 

with 6* : M \ {0} X M ^ M X M and e{x, y) = (x, |) such that 

With this adjustment the last step of the shearlet transform can be obtained with a standard 
inverse fast Fourier transform (similar as in our implementation). Unfortunately this is no 
longer related to translations of the shearlets in the time domain. 

3. Having a look at the Preprint of this paper, D. Labate made us aware of his Preprint \23^ 
where he and K. Guo proposed an interesting, new construction of smooth diagonal shearlets in 
contrast to our continuous diagonal shearlets ip^^^ . We implemented both constructions and 
as the smooth approach is promising for theoretical purposes the differences in our applications 
were negligible. 

3 Shearlet Regular izat ion in Image Segmentation 

In this section, we incorporate shearlets into a convex multi-label segmentation model. We 
start by explaining the convex relaxation model in matrix-vector notation. The approach is 
based on [SI IMlllS] and is closely related tol2l[l6l|26l[2Zl[39lllIl[lZ!- 

Segmentation Model. We segment an image by assigning different labels to different areas 
of the image based on the gray or RGB value of the pixels. We assume that the respective 
values are given in a codebook. In the following, we restrict ourselves to gray value images 
and will comment the slight adjustments for RGB images later. For a simple matrix-vector 
notation we assume the image to be column- wise reshaped as a vector. Thus, having an image 
/ G M^^^ we obtain the vectorized image vec(/) G where we retain the notation / for 
both the original and the reshaped image. Let g G N be the number of clusters which will be 
labeled by labels in C := {1, . . . , q}. 
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For a given image / G and a codebook c = (ci, . . . , Cg) we want to find a labeling vector 
n* G {0, that provides a "good" segmentation of /. For the interpretation it is more 
convenient to take u* G {0, i^^^^^i as q layers of an image in {0, 1}^^-'^. Usually for each 
pixel only the entry in one layer is 1 and the entries in the other layers are 0. Consequently 
the pixel gets the label according to the index of the respective layer. As a relaxation we 
allow non-zero entries in every layer but restrict the sum over all layers for each pixel to be 
equal to 1, i.e., u* E [0, 1]'^^" and Efc=o«*[^ + ^^^] = 1 for all 1 < r < N"^. The respective 
label is chosen according to the index of the largest entry for each pixel. 
The following variational approach finds u* as the minimizer of the functional 

F{u) = {u, s) + X^{u) subject to u £ C 

where 

q-l 

C ■.= {xe [0, 1]''^' : ^ x[r + kN^] = 1 for ah 1 < r < N^}. 

k=0 

The constraint guarantees that the sum over all layers for each pixel is equal to 1. We will 
refer to the first term as the data term and to the second as regularizer and A > as the 
regularization parameter. 

The relation to the given data / is contained in s G M'^^ which penalizes a certain distance 
between / and the codebook c. In each layer i of s we subtract the codebook value q from 
/ - in formulas: 

s[r +{k- l)N'^] = \\f[r] - c[k]\\P for 1 < r < iV^ I < k < q. (6) 

Thus, for big values in s the data term gets small for small values in u. On the other hand for 
small values in s we can chose values in u close to 1 what is needed for the sum to be equal 
to 1. The trivial solution u = is excluded due to the constraint u £ C. The regularization 
term is chosen to provide "smooth" layers. 

With the regularization parameter one can weight between data exactness and smoothness. A 
meanwhile frequent choice for the regularizer is the discrete TV-functional in its isotropic or 
anisotropic form. Note that for the anisotropic choice graph cut algorithms can be efficiently 
applied, see, e.g., [HP and [36] for a convex relaxation approach through continuous max- 
flow. This is also true for penalizers involving anisotropic variants of NL-means. While 
the TV-functional is well-suited for cartoon-like parts, NL-means are a powerful tool for 
restoring textured regions in images [6l [20l [T9] . We will use the TV- functional and NL-means 
regularizers in an isotropic form for numerical comparisons in Section [4j 
In this paper, we propose to minimize the functional 

F{u) = {u, s) + \\A{Ig S)u\\i + Lc{u), 

where Iq is the identity matrix in M'^'^'J and S is the shear let-transform applied to a vector in 
M'^^ in the sense of ([s]). In particular we have that 

S-S = I^2. (7) 

Finally, Iq S simply applies this shearlet-transform to each "layer" of u separately. The 
indicator function lq is defined by 



I'ciu) 



0, n G C 
oo, u ^ C. 
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Observe that our "regularization parameter" A is a diagonal matrix since we want to apply 
different parameters to different scales of the shearlet transform. With these notations we can 
choose different parameters for each scale and shear and even for each layer. In the following 
we will use the same parameter for each shear but different parameters for each scale. Thus, 
we simply write A = (Aq, Ai, . . . , Xjg) where Aq is the parameter for the "low frequency" and 
Ai, . . . , Ajo are the parameters for the different scales. 

Summing up, finding an appropriate labeling vector u* for the image / is equivalent to solving 
the problem 

u* = argmin{('u,s) + ||A(/g (g) S)u\\i + iciu)}. (8) 

Alternating Direction Method of Multipliers. To solve the problem in ([s]) we want to 
use the ADMM which applies in general to 

argmin {Gi{x) + G2{y)} subject to Ax = y 



with proper, closed, convex functions d -.W^ ^ MU{oo}, G2 : ^ MU{oo} and A £ M"^"'. 
The ADMM splits the problem into the following smaller subproblems which have to be solved 
alternatingly. 

General ADMM Algorithm: Initialization: e and 7 > 0. 

For i = 0, . . . repeat until a convergence criterion is reached: 

^(i+i) ^ argmin jGi(x) + — + Ax - y'^'^Wl] , 

y{^+l) = argmin |G2(y) + + Ax^-^^'^ - yfA , 

6(^+1) = 6W+^^(m)_y(m)_ 

To apply this general ADMM algorithm we rewrite the problem ([s]) as 

argmin {{u, s) + A||u||i + lc{w)} subject to {Iq (g) S)u = v, IqN2U = w (9) 

where r] is the number of scales and shears in the shearlet transform. Using Gi{u) := {u, s), 
G2{v,w) := A||u||i + ic{w) and A := f'^i '^'^\ -^g obtain the following algorithm: 

V ^qN2 J 

ADMM Algorithm for Initialization: g wniN^^^io) g r^a^'^^W ^ ]^{ri+i)qN^ and 
7 > 0. 

For i = 0, . . . repeat until a convergence criterion is reached: 

^(.+1) = argmin I (n,.) + + An - (""t) g] , (10) 

(jm)) = argmin (aIIHIi + ic{w) + + Au^^-^'^ " nO - (H) 
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The last step (12) is a simple update of b and can be split into two parts: one for "w" and 
one for "ti;". We obtain 



We have to comment the first two steps. Due to the differentiability of ( 10 ) its solution follows 
by setting the gradient of the functional on the right-hand side to zero. Thus, the minimizer 
is given by the solution of 



= s + {au + 6« 



A^Au = A 



— 7S 



By ([T]) we obtain 



A^A = {I^C^S^,IgN^] 



Ig(g)SjS+I, 



qN^ — "^IqN^ ) 



so that no linear system has to be solved and we get simply 



2 



(i) 



W 



-7s 



1 r 



which means that we have just to apply an inverse shear let transform. 



The second step (11) can be minimized for v and w separately. Thus, we have 

argmin|||At;||i + -^||^; - {b^^ + {Iq® S)u^'+^^) 



argmm 7 



1 + -^Wv- Qvh 



■-■gv 



and 



w 



(i+l) 



S'yA{gv) (componentwise soft-shrinkage) 



■gmin{.c(^) + ^\\w- (6« + V2n(^+i)) \\l} 



arg 



= Pc{gw) ("pixelwise" orthogonal projection onto C). 
The soft- shrinkage S\ with a threshold A is defined as 



Sx{x) := { 



X — X for 2; > A, 

for 2; G [-A, A], 

X + A for 2; < —A. 
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The set C is a subset of the hyperplane 

{x e R« : ^ Xfc = (x, Ig) = 1} C M'' 
k=l 

The orthogonal projection of g gM'^ onto the hyperplane is given by 

This can easily be constructed geometrically by calculating the intersection point of the 
straight line through g with slope Ic and the hyperplane. It cannot be guaranteed that g lies 
in C, i.e., some components could be less than zero. The projection onto C can be achieved 
by setting these components to zero and project again on a lower-dimensional space and so on. 

Finally we remark that primal-dual methods as in [10^ HH] could be applied instead of the 
ADMM. However, due to the Parseval frame property which preserves us from solving a linear 
system of equations in the first ADMM step, the primal-dual methods are not superior here. 

4 Numerical Examples 

In this section we demonstrate the performance of our algorithm by numerical examples. 
We will compare our method with minimization models using the same data term, but a 
different regularizer namely the classical TV regularizer and an NL regularizer. While the 
TV-functional is well-suited for cartoon-like parts, NL-means are a powerful tool for restoring 
textured regions in images [U EOj [21] . 

There exist different possibilities to apply the TV-like functional in segmentation. Here we 
focus on the method proposed in |32j which can be written in matrix- vector notation as 

N N ^ 

TV(n) := II |(/, C3 V)u\ ||i = J] |(/q ® V)u\{j) = ^ {\Vui{j)\^ + ■■■ + |Vtx,(j)|') ' , 

j=i i=i 

see [12]. As mentioned one could use alternatively + " ' " + |V^q(j)l)) where we 

have not seen a visually different result in our applications. 

Concerning the NL regularizer, we have used the approach in \19\ HI] with respect to the q 
labels. This means that we replace the discrete gradient matrix V by a matrix D which is 
constructed as follows: 

Initially, we start with a zero weight matrix w. For every image pixel i = (11,12), we compute 
for all j = (ii, j2) within a search window of size w x zu around i the distances 

da{i,j)-= Yl Yl Ga(^l,^2)(/(il + ^l,^2 + ^2)-/0l + ^l,J2 + ^2))^ 

where Ga denotes the discretized, normalized Gaussian with standard deviation a. We refer to 
p as the patch size. Then, for given m we select k < m so-called 'neighbors' j ^ i ofi for which 
da{i,j) takes the smallest values and set w{i,j) = w{j,i) := 1. By setting w{i,j) = w{j,i) it 
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happens that several weights w{j, •) are aheady non-zero before we reach pixel j. To avoid 
that the number of non-zero weights becomes too large, we choose only k = min{m, 2m — 1} 
neighbors for I being the number of non-zero weights w{j, •) before the selection. Now we 
construct the matrix D G W^^'^ with m = 2m so that D consists of m blocks of size N x N, 
each having —1 as diagonal elements plus one additional nonzero value 1 in each row whose 
position is determined by the nonzero weights w{i,j) and maybe some zero rows. The fol- 
lowing parameters were used throughout the experiments: p = 5,w = 15,a = 2,m = 5. 



First we present an artificial example. In Fig. [5] an original — 1 grid image and its noisy 
version with white Gaussian noise of standard deviation 0.2 are shown. The noisy image was 
segmented by minimizing the proposed convex functional with TV regularization (bottom 
left) and shearlet regularization (bottom right). As codebook we have used the exact values 
and 1. While the first method fails to reconstruct the grid, our shearlet segmentation method 
works perfect. An exact reconstruction can be also obtained by applying the NL- means 
regularizer. However in its isotropic version the later method requires more time than the 
shearlet approach. 

Next we want to segment the cartoon type image in Fig. [6| The top left image (Fig. |6(a)| ) 
shows the original cartoon image with four constant regions. The image was taken from [7J. 
In the top right image white Gaussian noise of standard deviation 0.1 was added (Fig. [6(b)] ). 
The bottom row shows the difference between the segmented image and the original one. For 
the left image (Fig. 6(c) ) the TV regularizer was used and for the right image (Fig. |6(d)[ ) 



the shearlet regularizer. In both images a white pixel represents the correctly labeled pixels 
and the black pixels indicate false labeled pixels. Both methods have problems with correctly 
labeling the borders of the circles and the TV method also fails on assigning the correct labels 
for the diagonal border of the triangle whereas the shearlet regularizer labels all these points 
in the right way. 

To illustrate the good performance of the shearlet regularizer for curved structures com- 
pared to the TV regularizer we provide another example using the april calender sheet of 
www.mathe-kalender.de, i.e., a snippet from the right bottom corner (see Fig. |7(a)| ). For 
the segmentation we added again white Gaussian noise of standard deviation a = 0.1 in Fig. 
|7(b)[ The six different reference colors for the codebook were chose manually. In Fig. |7(c 



we show the result using the TV regularizer and in Fig. |7(d)| we used the shearlet regularizer. 
Our shearlet method preserves the stripes better than the first method. Adding more noise 
the results using TV become worse whereas the results using shear lets remain similar. 



Last we want to segment the real image in Fig. [8| where 8(a) is the original RGB image of 
a clown-surgeon fish and in |8(b)| white Gaussian noise of standard deviation 0.2 was added. 
As codebook we use the matrix 

^0.7451 0.1843 0.3686 0.8353^ 
0.8314 0.2784 0.5569 0.7333 
^0.8196 0.2275 0.6353 0.3020^ 

where each row represents a color channel and each column stands for one label. The colors 
were chosen manually according to the different regions of the image. Compared to the 
segmentation model above for gray valued images only the computation of s in ([6]) has to be 
slightly modified. Assuming that every pixel in the image / consists of a three element vector 
(/ij) /G) Jb) we compute s as the norm of the difference between the vectors {/r, fc, fs) and 
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(a) Original image 



(b) Noisy image 




Figure 5: Segmentation of a noisy grid image with TV regularization and shearlet regulariza- 
tion. 
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(a) Original cartoon image 



(b) Cartoon image, noise a = 10% 



(c) Difference between original and seg- 
mented image using TV regularized seg- 
mentation, 7 = 1, A = |, 100 iterations 



(d) Difference between original and seg- 
mented image using shearlet regular- 
ized segmentation, 7 = 1, A = 
^(0,0.1,0.2,2.20), 50 iterations 



Figure 6: Four-class of a cartoon image 
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(a) Original image with colored stripes (b) Image with colored stripes, noise a = 

10% 




(c) Segmented image using TV regularizer, (d) Segmented image using shearlet regular- 

7 = 1, A = ^, 100 iterations izer, 7 = 1, A = ^(0,0.1,0.2,0.4), 50 

iterations 



Figure 7: Segmentation of stripes 
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(chjCcCb) - in formulas: 



s[r + {k-l)N^ 





for 1 < r < iV^ 1 < A; < 





(a) Original image 



(b) Noisy image 



Figure 8: Original RGB image of a clown-surgeon fish and its noisy version with white 
Gaussian noise of standard deviation 0.2. 

In Fig. [9] we compare the segmentation results for both images when applying different 
regularizers in the minimization functional. 



The first row (Fig. |9(a) and |9(b) ) was segmented using the TV regularizer. The second row 



(Fig. 9(c) and 9(d)[ ) is "obtained with the NL regularizer. The last row (Fig. [9(e)] and [9(f)t 
shows our shearlet regularization. The parameter were chosen such that we get visually the 
best result. For the original image the results are rather similar with slight difficulties for 
the TV-model to segment the curves in the front of the fish. For the noisy image the results 
for the TV-based method get worse. As expected for different parameters either noise and 
structure is preserved or the noise vanishes in the segmentation process but we loose some 
structure, too. However, the TV method is the fastest one. The NL-means regularizer does 
a good job, but requires more computational time than our shearlet method and tends to 
preserve tiny structures. Our shearlet method segments this kind of curved structure very 
well. 

In summary we have seen that shearlets can be successfully applied for the segmentation of 
curved textures in conjunction with a convex multi-label model and the ADMM algorithm. 
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